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We have studied the contact network properties of two and three dimensional polydisperse, 
frictionless sphere packings at the random closed packing density through simulations. We 
observe universal correlations between particle size and contact number that are independent 
of the polydispersity of the packing. This allows us to formulate a mean field version of the 
granocentric model to predict the contact number distribution P{z). We find the predictions 
to be in good agreement with a wide range of discrete and continuous size distributions. 
The values of the two parameters that appear in the model are also independent of the 
polydispersity of the packing. Finally we look at the nearest neighbour spatial correlations 
to investigate the validity of the granocentric approach. We find that both particle size and 
contact number are anti-correlated which contrasts with the assumptions of the granocentric 
model. Despite this shortcoming, the correlations are sufficiently weak which explains the 
good approximation of P{z) obtained from the model. 
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1. Introduction 

The question of how spheres pack together has been of interest to scientists for 
centuries In the context of amorphous materials the jamming transition of ran- 
dom close packed spheres are of particular interest and have been a substantial area 
of study in recent yearsj3-0]- At the jamming point, which corresponds to a critical 
packing density, the packing makes a sharp transition towards a mechanical stable 
state. At this point the isostatic condition requires that the average coordination 
number of the packing is twice the number of dimensions of the system Q]. This 
density is referred to as the random close packing density (pRCP- 

The packing of equal sized spheres in disordered configurations have a long 
history [8| and the value at which (/>_rcp is reached for these packing is well stud- 
ied in both experiment [1, 0] and simulations 0, U], though 4>rcp has been 



shown to be dependent upon the history of the packing and the packing protocol 
used0, More prevalent in nature, though not as a widely studied, are packings 
with a distribution of sizes. Experiments and simulations of binary mixtures|13l- 



16[ and continuous size distributions'l7H22l] have investigated the value of (j)RCP 



and found that it increases with polydispersity. There have been some simula- 
tion and experimental studies on the contact properties of polydisperse, disordered 



packings|18l. I23l . |2j] and recently the granocentric model has been proposed to 



predict the local packing structure at (/'ijcp |21|, l2j, l25[ in three dimensions. 
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Figure 1. Visualisation of soft sphere packing simulation at ipRCP with a lognormal distribution of radii. 
The spheres are coloured with a spectrum going from blue to red according to size with blue correspond 
to the smallest particles. 

In this report we investigate the correlations between size and contact number 
of particles in polydisperse packings at 4>rcp in two and three dimensions for a 
wide range of size distributions. Our key finding is the existence of universal corre- 
lations between size and contact number that is independent of the polydispersity. 
This empirical result allows us to formulate a mean field approach based on the 
granocentric model that yields excellent agreement with our data. One of the key 
assumption in the granocentric model is the lack of spatial correlations of both 
size and contact number of the particles. Our measurements of nearest neighbour 
correlations show that this assumption is violated. In general, the average contact 
number and the average size of neighbouring particles do not correspond to the 
global mean of contact number and size. In 3D packings larger particles are sur- 
rounded by smaller particles and vice versa. Moreover, particles with few contacts 
are neighbouring particles with many contacts. Nevertheless, these correlations are 
weak enough so that the predictions we obtain from the granocentric model agree 
well with our data. 



2. Simulations 

We model the disordered packings at (pRCP through simulation of soft spheres. 
These are frictionless spheres that interact through purely repulsive body centred 
forces, which can be written as a function of the overlap between two particles in 
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Figure 2. Continuous size distributions used to create the soft sphere packings. The distributions are: (Q) 
lognormal aA = 0.40; ( ) Gaussian cta = 0.44; (□) uniform cta = 0.44. The open symbols represent the 
original size distribution and the closed symbols represent the size distribution once rattlers are removed. 

contact. The overlap is 



where Ri and Rj are the radii of spheres i and j and dij is the distance between 
the respective centres of the spheres. The interaction potential of the spheres is 



These interactions are harmonic with a spring constant k. The spheres have their 
radius drawn from a set size distribution and are placed at random in a three 
dimensional periodic cell. The radii of the spheres are then rescaled such that 
the desired packing fraction (/> is reached. A conjugate gradient method is then 
used to minimise the overlap between spheres and hence the the total energy of 
the packing [2^. The simulation is halted when the packing is in a local energy 
minimum and in mechanical equilibrium. 

In general, the isostatic condition can be shown by the following argument. If 
there are D dimensions with soft particles, at (f)RCP there will be an average 
number of contacts {z) (where (.) denotes an average over all particles). Therefore 
in total there will be N{z)/2 contacts in the packing since every contact is shared 
by two particles. For mechanical stability all the contact forces need to balance on 
each particle [3], which leads to matching the ND degrees of freedom with the the 
contact forces giving 





ND=m 
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Figure 3. (a) ct>RCP versus the standard deviation of the surface area distribution a a for a variety of 
size distributions in three dimensions: (v) monodisperse; (0) bidisperse; (□) uniform; ( ) Gaussian; (Q) 
lognormal. ipRCP including all particles is plotted with open symbols while solid symbols correspond to 
<t>RCP with rattlers omitted. The dashed line indicates the average 4>rcp with rattlers omitted for all size 
distributions. Inset: Percentage of rattlers at <j>RCP in three dimensions versus the standard deviation of 
the surface area distribution tr^. (b) ipRCP versus the standard deviation of the radius distribution in 
two dimensions with (open symbols) and without rattlers (closed symbols). The dashed line is the average 
4>RCP with rattlers omitted for all size distributions with an > 0.1. Inset: Percentage of rattlers at (pRCP 
in two dimensions and the standard deviation of the radius distribution ctr. The symbols correspond to 
the same data as in Figure [Sfa). 



The isostatic point Zc, equivalent to (j)RCP, which is defined as when 



in 3 dimensions, 
in 2 dimensions. 



(4) 



from Equation While globally these mechanically jammed states are con- 
strained to have {z) = Zc, there is a distribution of contact numbers for particles. 

In general particles that have less than D + 1 contacts cannot be locally mechan- 
ical stable. For 3D that means all particles with less than 4 contacts and 2D, all 
particles with less than 3 contacts are locally unstable. These particles are called 
rattlers and their contribution to the contact number analysis is omitted as their 
contact number is ill-defined. In a recent publication [27] it was shown that pack- 
ings which satisfy Equation ^ can still be unstable under shear, though this effect 
is only pronounced for systems with a low number of particles in the packing. This 
effect is negligible in our simulations. 

Each of the simulated packings has 16384 particles with various different size 
distributions, with up to 500 realisations in three dimensions and 50 realisations 
in two dimensions for each size distribution. An example of a sphere packing is 
shown in Figure [TJ A variety of size distributions are created including discrete 
size distributions of monodisperse and bidisperse spheres, where there is a 50-50 
mixture with a size ratio 1:1.4, and continuous radius distributions such as the 
lognormal distribution, Gaussian distribution and uniform distribution, which are 
plotted in Figure [2j Packings at 4>RCP found by starting with a packing density 
above (pRCP which is lowered until the average contact number (z) is within the 
range Zc + 0.05 > (z) > z^. 
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Figure 4. Variance of the contact number distribution a\ versus the variance of the area distribution cr^. 
The dashed line corresponds to a linear fit to the data: (t^ = 1-60 + 8.09(7^. The closed symbols are the 
predictions from Equation 1301 Plotted in the inset is the variance of cr^ versus the radius distribution (t|j 
for two dimensional packings. The dashed line corresponds to a linear fit to the data: cr^ = 0-61 + 3.52(t|j. 
The closed symbols are the predictions from Equation 1331 The data is labeled as in Figure [3ja) . 

2.1. Properties of Polydisperse Packings 

Next we define apt as the normalised standard deviation of the P{R) distribution 
where rattlers have been removed, 



It was also found useful for packings in three dimensions to define the standard 
deviation of the corresponding normalised surface area distribution as 

The size distribution affects the packing density at which the isostatic point is 
reached jl7]. As the width of the size distribution is increased, (pRCP becomes larger 
because smaller particles are able to fit between the interstices of larger particles 
in contact as seen in Figure [3]^a) for three dimensions and in Figure [3|^b) for two 
dimensions. This also results in an increase of rattlers [l9^, as shown in the insets of 
Figures [3)[a) and (b). (j)RCP only depends strongly on the width but not the shape 
of the size distribution. As the size distribution becomes wider the percentage of 
rattlers increases. For polydispersities with a large population of small particles 
such as the uniform distribution this results in an increase of rattlers of up to 50%, 
though (pRCP is only slightly affected. 

Also plotted in Figure [3]J^a) and Figure [3^b) is (j)RCP when the volume of rattlers 
is excluded. This (pRCP with rattlers omitted is found to be a constant that is 
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independent of the size distribution in three dimensions, where the average 4'rc'p 
with rattlers omitted is 0.621 it 0.003. In two dimensions for aji > 0.1 the (pRCP 
with rattlers omitted is also constant and independent of polydispersity with the 
average (pRCP with rattlers omitted equal to 0.803 it 0.002. Two dimensional disc 
packings with aR < 0.1 partially crystallise [28] which leads to the increase in in 
(pRCP and rattlers for packings in two dimensions with u/j < 0.1. 

Changing the polydispersity also affects the contact properties. As shown in 
Figure [H changing the width of the size distribution affects the variance of the 
contact number distribution cr^- The standard deviation of the contact number 
distribution az is defined as az = \J {z^) — {z)"^- Broader size distributions results 
in broader contact number distributions. This trend is independent of the type 
of size distribution in both two and three dimensions. We find that o"| increases 
linearly with aj^ in two dimensions and a\ in three dimensions. For two dimensional 
cellular structures a corresponding relationship between the standard deviation of 
the size distribution and the standard deviation of the number of cell faces has 
been observed 2^31]. While the width of the contact number distribution is set 
by the width of the size distribution only, its shape does depend on the particular 
size distribution as can be seen in Figure [9l 



3. Local correlations in polydisperse packings 

3.1. Contact number and size correlations in three dimensions 

While a large body of literature on random packings is devoted to the bulk prop- 
erties of mono- and bi-disperse packings near the jamming transition [3, 0, 0, [13], 
important results on the local structure in polydisperse packings have emerged 
only in recent years . The pioneering work by Clusel et al. established a 



link between the size distribution and the local structure of the packing. We have 
expanded on this work by investigating how the correlations between particle size 
and contact number depend upon the polydispersity of the packing. 
The average contact number for particles of a given size is defined as 



(z|x) = J]zP(z|x), (7) 



where P{z\x) is the contact number distribution for particles of a given size x. 
The average contact number for particles of given size x at (pRCP is plotted for 
a wide range of size distributions of different widths and shape. We scaled the 
data in three different ways, in terms of the normalised radius, normalised surface 
area and normalised volume in Figure [5^a), (b) and (c) respectively. In the three 
scalings {z\x) for all size distributions and polydispersities follow similar trends. 
Namely, larger particles have more contacts on average. This can be explained in 



the context of the granocentric model [2l|, [24] which stipulates that larger particles 



have more solid angle available to accommodate neighbouring spheres. Surprisingly, 
these correlations are independent of polydispersity. The best collapse is observed 
when the scaling is in terms of the normalised area 
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Figure 5. The average of the contact number distribution for a particle of a given size for the six different 
size distributions at <f>iiCP'- iy) monodisperse; (0) bidisperse, radius ratio 1:1.4; (□) uniform cta = 0.23; 
(a) Gaussian cfa = 0.27; (O) lognormal a a = 0.40; (<) lognormal oa = 0.61; (*) lognormal aA = 0.72. 
We present three different scahngs: (a) in terms of the normaUsed radius r; (b) in terms of the normalised 
area a; (c) in terms of the normalised volume v. The data are plotted over a range that illustrates the 
quality of the collapse. 



as shown in Figure [5jb). This cohapse of the data is well described by a linear fit 

(z|a) = (z)+7(a-l), (9) 

which is plotted in Figure [6l The form of Equation Q ensures that the isostatic 
constraint 

oo 

{z) = j {z\a)P{a)da = 6, (10) 



is satisfied. The fitting parameter is found to be 7 = 3.032 it 0.004. The contact 
number average {z\a) for the discrete distributions (monodisperse, bidisperse) has 
the same value as that of a particle of the same size in the continuous distributions. 
Figure [5]^b) shows that at (pRCP, the relationship between z and a is universal and 
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Figure 6. The average contact number for particles of a given area a at <f>RCP in three dimensions for 
all (+) monodisperse, (0) bidisperse, (□) uniform, ( ) Gaussian and (Q) lognormal size distributions at 
all the widths cr^ we have considered (see Figure |4]l. The solid red line is a fit to Equation l(9]|- Inset: the 
average contact number for particles of a given radius r in tv^o dimensions at 4'RCP- The solid red line is 
a fit to Equation II22I I. 



independent of size distributions. This suggests that the local contact properties of 
a particle only depends upon its surface area. This result is similar to that observed 
in two dimensional disc packings 33|], which will be discussed in Section [3.31 

It must be noted that the average {z\a) > 4 since we omit rattlers from our 
analysis. Also, for large values of a the scatter in {z\a) is much larger due to lower 
statistics. For all equations fitted and figures plotted with the exception of Figures 
El [71 m [lOl [I2l [I3l data binned with less than 100 particles are omitted. 

In Figure[71 a number of different contact number distributions P{z\a) are plotted 
for given intervals of a. This figure demonstrates that the P(z\a) distributions are 
independent of shape and width of the size distribution. This confirms what is 
suggested in Figure [6] - namely that the contact number distribution for a particle 
in a packing at (pRCP does not depend on the global size distribution of the packing 
but on the size of the particle only. 



3.2. Mean field granocentric model in three dimensions 

We have shown the contact number distributions P(z\a) for particles of a given 
size do not depend on the global size distribution of the packing but only on the 
size of the particle in question. This result allows us to formulate a mean field 
granocentric model that is similar in spirit to the one by Newhall et al. [3J] who 
investigated size-topology relations in tessellated packings. 

Here we use a mean field approach that allows us to predict the correlations 
between size and contact number. In contrast to the original granocentric model 



2l[ | we explicitly exclude rattlers since their contact number is ill-defined. 

Since local correlations are independent of the size distributions, we consider a 
particle of a given radius Rc and then make a mean field assumption that all the 
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Figure 7. The contact number distribution for particles of a given size P{z\a) for four size distributions 
at rpRCP in three dimensions: (v) monodisperse; (□) uniform aA = 0.44; ( ) Gaussian ua = 0.44; 
(O) lognormal cta = 0.40. Plotted here is a selection of the P(z\a) for 6 different intervals a with the 
P(z\a) shifted for clarity. Plotted in order of lowest to highest is 0.475 < a < 0.525; 0.975 < a < 1.025; 
1.475 < a < 1.525; 1.975 < a < 2.025; 2.475 < a < 2.525; 2.975 < a < 3.025. The solid red line is the 
model prediction of P{z\a) from Equation I I16I I. 



particles surrounding it are of average radius {R). If this particle of size Rc is in 
contact with another particle {R) it will subtend a solid angle of the central 
particle, which is given by 



n{R,, {R)) = 2n{l- ^ + 1^ I • (^1) 

Since we have shown that the proper scaling of correlations between size and con- 
tact number is in terms of a we rewrite Equation (llip accordingly with all contact- 
ing particles now assumed to have an average radius ^/{R?)\ 

J)(a) = 2^(l-— )/^Jl + ^l . (12) 
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Having obtained a value of the solid angle subtended by a single contact, the 
maximum number of contacts is simply 

A correction must be made to Zmax to account for the interstices, similar to the 
familiar sphere kissing problem for monodisperse spheres where only 12 spheres 
can be in contact with a central sphere even though there is sufficient solid angle 
to fit 14 spheres [IH]. 

A prefactor a is introduced into the model to limit the maximum number of 
contacts: 

2a 

Zmax{a) = p , ^ • (14) 



In order to recover the known result of the kissing problem for monodisperse 
spheres, the value of a would have to be 0.8708. In our model, however, the value 
of a will turn out to be less than that due to additional constraints. 

Following the granocentric approach [21l], we now make an ansatz that the dis- 
tribution of the number of particles in contact with a particle of size o is given by 
a binomial distribution. 

P{z\a) = B{z; Zmax{a),p), (15) 



where B{z] Zmax{o),p) is a binomial distribution with the maximum number of 
trials, that is the number of times in which a particle can attempted to be placed 
in contact with the particle of size a, given by Zmax and p is the acceptance 
probability that a particle will be in contact. The probability p is the other free 
parameter in the model. 

In order to omit rattlers we truncate the binomial distribution for z < 4 by 
including a Heaviside function H{z — 4) and a normalisation constant C, so P{z\a) 
becomes, 

Piz\a) = B'iz;Z^axia),p), (16) 
= C .^ p^{l-pf--^H{z - 4). 

Note that this is in contrast to the original granocentric model which did not 
exclude rattlers. This allows us to make a prediction for the contact number average 
for a given particle size, 

{z\a)= zPiz\a), (17) 

z=A 
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Figure 8. (a) The average of the contact number distribution for particle of a given size for five different 
size distributions at ipRCP- The model prediction from Equation l|17| l is plotted as the solid red line, (b) 
The variance of the contact number distribution for particle of a given a. The model prediction from 
Equation l|18| l is plotted as the solid red line, (c) The ratio of the variance to the average of the contact 
number distribution for particles of a given size a. The dashed line denotes the value of the acceptance 
probability p as found from the data. The size distributions and symbols plotted in all three panels are the 
same as in Figure[7]with the addition of: (0) bidisperse, radius ratio 1:1.4 and (*) lognormal (Ta = 0.72. 



and the corresponding variance of the contact number for a given particle size 

(4l«)= E (18) 

z=4 

as weh as a prediction of the global contact number distribution of the packing 

oo 

P{z) = J P{z\a)P{a)da. (19) 



Given a size distribution P{a) and using Equation ()16p . a prediction for the 
contact number distribution can be made for any packing at (jjRCP- The acceptance 
probability p can be determined through a property of the binomial distribution. 
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If X ~ B(n,p) is a random variable from a binomial distribution B with n trials 
then the mean is given by 



E[X] = np, 



and the variance is given by 



yar[X] = np{l — p). 

Therefore the ratio of the variance to the mean of a binomial distribution is a 
constant given in terms of p, which in the context of our model is given by 



{z\a) 



Equation (jl6p is a truncated binomial distribution and therefore the ratio '^jjj^"' 
plotted in Figure [8{c) is only expected to reach a constant at sufficiently large 
values of a where the truncation becomes negligible. Indeed, for a > 2 the ratio 
plateaus at 0.342 ± 0.006 which corresponds to p = 0.658 it 0.006. 

After obtaining the probability p directly from the data we can fix the second 
parameter a by imposing the constraint as stated in Equation (llOp . namely that the 
global average contact number of the packing (z) must be equal to 6. This results 
in a taking a value of 0.625. Surprisingly, a does not depend on polydispersity, 
therefore the two free parameters of the model, a and p, can be fixed for all size 
distributions at 4>rcp- This may be related to the fact (pRCP without the rattlers, 
which are explicitly omitted in this model, is a constant (Figure [3]) . 

The constancy of the two parameters is a significant simplification to the original 
granocentric model, where the acceptance probability and maximum solid angle 
need to be determined for each polydispersity separately. 

Comparing the prediction of the average {z\a) from Equation (I17p with the data 
as shown in Figure El^ a) , we see that the model is in good agreement with the data 
over a large range of a. Only for large values of a it deviates slightly. Similarly, 
the model prediction of the variance {cr'^\a) (Equation (llSp ) agrees well with the 
data as shown in Figure El^b) . Note that the staircase structure exhibited by the 
model in Figure [8]^ a) -(b) is due to the discrete nature of the binomial distribution. 
Finally, we can compare the prediction for the contact number distribution P{z) 
from Equation (|19p with the data for a wide range of size distributions as shown 
in Figure [9l As mentioned earlier, the same parameters a and p are used for all 
polydispersities. 

For the continuous size distributions the agreement is excellent, with slight de- 
viations in the tails, while for the discrete size distributions the model fails to 
reproduce the tails of the distribution. This is a consequence of the a parameter 
which takes on a value that limits the number of contacts to less than the maximum 
number of contacts allowed by geometry. 

For example, to recover the maximum number of contacts in monodisperse pack- 
ings, a would need to be 0.8708. This would violate the constraint imposed by 
Equation (|10p as the probability p is set by the data and cannot be adjusted. As 
the size distribution becomes wider this discrepancy becomes less pronounced. 
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Figure 9. The contact number distribution for five different size distributions at <j>RCP- The P(z) are 
shifted for clarity. The open symbols represent the simulation data and the closed symbols represent the 
model prediction from Equation I I19I I . The size distributions and symbols plotted are the same as in Figure 



3.3. Mean field granocentric model in two dimensions 

We now use the same approach for 2D polydisperse packings. In Figure [TOT a). 
(b) and (c) the average contact number for particles of a given size x at (pncp 
is plotted and scaled in terms of the normalised radius, normalised surface area 
and normalised volume. In the three scalings the {z\x) for all size distributions 
and polydispersities follow similar trends, namely that larger particles have more 
contacts on average. Similar to the three dimensional case, the best collapse of the 
data is found when the scaling is 

^ (21) 



as shown in Figure fTOTa). Therefore, the proper scaling variable for size-contact 
number correlations in 2D packings is x = r. The inset of Figure [6] exhibits a 
similar collapse of the average contact number for particles of a given radius {z\r) 
for a wide range of size distributions to that found for {z\a) in three dimensions. 
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Figure 10. The average of the contact number distribution for a particle of a given size for the six different 
size distributions in two dimensions at <j>RCP- (0) bidisperse, radius ratio 1:1.4; (□) uniform an = 0.17; 
{<) lognormal cr/j = 0.10; ( ) Gaussian ct/j = 0.24; (v) lognormal ur = 0.35; (t>) lognormal an = 0.45. 
We present three different scalings: (a) in terms of the normalised radius r; (b) in terms of the normalised 
area a; (c) in terms of the normalised volume v. The data are plotted over a range that illustrates the 
quality of the collapse. 



Similar to Equation ([9]), the average contact number {z\r) is well fit by a linear 
function of the form 

{z\r) = {z)+-f2D{r-l), (22) 

with the fit parameter 72D = 2.023 it 0.007. 

In Figure [11] we plot the two dimensional equivalent of Figure [8l except with 
the size of the particle now represented by the normalised radius r instead of the 
normalised surface area a. The model prediction that appears in Figure [TTTa) for 
{z\r) is analogous to the model outlined in the previous section. The principle 
adjustment is that the maximum number of discs that can be placed in contact 
with a disc of radius r must now be expressed in terms of the available angle rather 
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Figure 11. (a) The average of the contact number distribution for particles of a given radius r in two 
dimensions for six different size distributions at <f>RCP'-(^) bidisperse, radius ratio 1:1.4; (□) uniform 
aji = 0.17; (a) Gaussian aji = 0.24; ( :) lognormal crji = 0.20; (t>) lognormal cr^ = 0.30; (O) lognormal 
an = 0.40. The model prediction is plotted as the solid red line, (b) The variance (cr^lr) for the same size 
distributions. The model prediction is plotted as the solid red line, (c) The ratio of the variance to the 
average of P(2|r). The dashed red line denotes the value of the acceptance probability p as found from the 
data. 



than the sohd angle. Equation (jl4p now reads as 

avr 

Zmax{r) = _ / N ) (23) 
[—r ) 

which affects the number of trials of the binomial distribution. Also rattlers, par- 
ticles with less than 3 contacts in 2D packings, are excluded from the binomial 
distribution. 

The acceptance probability p can be determined in the same fashion using the 
ratio between the variance of the contact number for particles of a given radius 
(o"^|r) and the corresponding average {z\r) as shown in Figure [TlTc). As before, 
the ratio appears to plateau for large particles which is consistent with a truncated 
binomial distribution. The value for p found from the data shown in Figure [TlTc) 
is 0.78 ± 0.02. 
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Figure 12. The contact number distribution for discs of a given radius r, P{z\r) for the same size distribu- 
tions as plotted in the the inset of Figure [TTI at t/iRCP- Plotted here is a selection of P{z\r) for five different 
intervals r with the P(z\r) shifted for clarity. Plotted in order of lowest to highest is 0.475 < a < 0.525; 
0.825 < a < 0.875; 0.975 < a < 1.025; 1.125 < a < 1.175; 1.475 < a < 1.525. The solid red line is the 
model prediction of P{z\r). 



Analogously to the 3D case, a is determined by the isostatic constraint which is 

oo 

(z) = J {z\r)P{r)dr = 4. (24) 


in 2D packings with the size distribution now given in terms of r rather than a. 
The value of a for two dimensions is found to be 0.894. The equivalent value a 
that would recover the correct answer for the kissing problem in two dimensions is 
1. 

The justification for this model is that the correlation between size and contact 
number is independent of of polydispersity, analogous to the results in three di- 
mensions. Similar to Figure [7] a number of difference size distributions are plotted 
for given intervals of r in Figure [T2l For each r interval plotted all the P{z\r) col- 
lapse independent of size distribution therefore validating the basis of the model. 
However, Figure [12] highlights some of the weaknesses of the model. For example, 



17 




z 



Figure 13. The contact number distribution for a number of different types of polydispersities at ipRCP 
in two dimension. The data from simulation is plotted as open symbols and the prediction from the model 
is plotted as closed symbols. The parameters p = 0.78 and a = 0.894 are used for all size distributions. 
The same size distributions are plotted with the same symbols as in Figure [TT] Distributions arc shifted 
for clarity. 



for the lowest interval of r where the data shows a range of contact numbers z 
but the model predicts that only 3 discs can fit around a disc of that size. This 
discrepancy at low r is due the contact limiting parameter a. For larger values of 
r the model prediction of P{z\r) is in better agreement with the data. 

The global contact distribution can then be predicted from the two dimensional 
equivalent of Equation (fT9]) , 



P{z) = / P{z\r)P{r)dr. 



(25) 



Figure [13] shows good agreement between the predictions and data for a wide range 
of size distributions. Similar to the results in three dimensions, the prediction is in 
closer agreement with the data for wider size distributions. 
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3.4. Size of a particle with contact number z 

3.4-1- Size of a particle with contact number z in three dimensions 

We have shown that {z\a), the average contact number for a particle of a given 
size is independent of the size distribution. However, it is important to emphasise 
that the converse is not true. In general, {a\z), the average area of particles that 
have z contacts, which is defined as, 

POO 

{a\z) = / aP{a\z)da, (26) 
Jo 

is not equal to {z\a). The bottom inset of Figure [T^ shows that {a\z) is not inde- 
pendent of size distribution. 

Continuous size distributions, lognormal and Gaussian are well approximated by 
a linear relationship 

{a\z) = 1 + Xiz- {z)). (27) 

This functional form ensures that '^z{a\z)P{z) = 1. When rescaled by the fitting 
parameter A, {a\z) collapses for all lognormal and Gaussian size distributions, as 
shown in the top inset of Figure [TH While the overall trend of {a\z) for lognormal 
and Gaussian distributions is linear there are deviations. 

Size distributions that lack tails, such as the discrete distributions, have a differ- 
ent functional form of {a\z) because of the large population of big spheres that can 
take a wide range of z as seen in Figure [71 This causes {a\z) to plateau at large z, 
as shown in the bottom inset of Figure [Hi In the limit of monodisperse packings, 
{a\z) = 1, which corresponds to A = 0. 

The linear relationship between size and contact number is similar to Lewis' law 
(sgI I for two dimensional cellular structures. 

The two different relationships between a and z arise from being calculated from 
two different conditional probabilities, the discrete distribution P{z\a) and the 
continuous distribution P{a\z), which are related by Bayes Theorem, 

P{z\a)=P{a\z)^. (28) 

Therefore, P(z\a) and P{a\z), in addition to being discrete and continuous distri- 
bution respectively, are related but not the same. Hence the fitting parameter A is 
not simply the inverse of 7. 

From Equation (p8]) . the two conditional averages {z\a) and {a\z) can be related 
by 

POO 

/ a(z|a)P(a)da = Vz(a|z)P(z). (29) 

Jo z=A 

Then substituting the linear fits of Equation ([9]) and Equation (j27p into Equation 
(j29p a relationship can be found between the width of the size distribution and the 
width of the contact number distribution, 



4 = TTY-A (30) 
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Figure 14. The average area of particles with a given contact number for four different size distributions in 
three dimensions at <j>jicp'- (<) lognormal cr^ = 0.10; (. ) Gaussian cta = 0.27; (O) lognormal cr^ = 0.40; 
(>) lognormal aj\ = 0.82. The solid lines are fits to Equation 1271 Top inset shoves the average area of 
particles with a given contact number for all (Q) lognormal and ( ) Gaussian size distributions at (pRCPy 
collapsed after rescaling by fitting the data to Equation I I27I1 . The dashed red line corresponds to a slope 
of 1. Inset of the top inset shows the fit parameter A as a function of a a- Bottom inset shows relationship 
between {a\z) versus z for the same size distributions as plotted in Figure |8] using the same symbols. 

While 7 is a constant, A clearly depends on the width of the size distribution aA as 
shown in the inset of Figure HM Although the functional form of A(o"yi) is unclear, 
we can substitute the values for A into Equation ()30p to compare with the data 
for cr^ versus a\ from Figure IH The agreement is good for broad distributions but 
less accurate for narrow ones where {a\z) is not well approximated by the linear fit 
(Equation [271). 

3.4-2. Size of a particle with contact number z in two dimensions 

Similar correlations are also observed in two dimensions. By defining the average 
radius for particles with a given contact number {r\z) as. 



we find that {r\z) is well approximated by a linear relationship similar to Equation 



The fits to this equation shown in Figure [15] agree well for the 2D Gaussian and 
lognormal disc distributions. Analogously to the 3D packings, {r\z) for size distri- 
butions without tails shows deviations from the linear fit at large z as shown in the 
bottom inset of Figure [T5l though it is less pronounced due to the smaller range 
of contact numbers in 2D packings. 

In the top inset of Figure [15] the correlations {r\z) for all lognormal and Gaussian 
size distributions are rescaled to highlight the agreement with Equation (j32p . The 




(31) 




(32) 
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Figure 15. The average radius of particles with a given contact number for four different size distributions 
in two dimensions at (pRCP- (*^) lognormal an = 0.10; ( ) Gaussian an = 0.24; (Q) lognormal aji = 0.30; 
(>) lognormal an = 0.45. The solid lines are fits to Equation 1321 Top inset shows the average area of 
particles with a given contact number for all (Q) lognormal and ( ) Gaussian size distributions at (pRCPy 
collapsed after rescaling by fitting the data to Equation I I32I1 . The dashed red line corresponds to a slope 
of 1. Inset of the top inset shows the fit parameter A as a function of a a- Bottom inset shows relationship 
between {a\z) versus z for the size distributions: {()) bidisperse; (□) uniform an = 0.23; (. .) Gaussian 
an = 0.27; (O) lognormal an = 0.35. 



fitting parameter, A2D is shown in the inset of Figure [15] as a function of the 
polydispersity ctr. 

Analogously to the 3D case, a relationship between the standard deviations of 
the size distribution and the contact number distribution can be written assuming 
that size-contact number correlations are approximately linear, 

2 72D ^2 /oo^ 

Comparison of this relation with the data is shown in the inset of Figure [H While 
this relation works well for broad distributions, narrow ones are not captured well 
for the same reasons as in the 3D case. Also, in 2D packings partial crystallization 
takes place for narrow size distributions. 



4. Nearest neighbour correlations in packings 

4.1. Nearest neighbour contact number correlations 

In the previous section we proposed a mean-field model based on the assumption 
that the packing is spatially uncorrelated. Specifically, we assume that the contact 
number and size of a particle is uncorrelated to the contact number and size of its 
contacting neigh ours. This assumption is implicitly made in many recent models 
that predict the density[ll|], contact distribution^, 34| and force networks [37i] in 
disordered packings. 
It is therefore worthwhile to investigate this in more detail. To do so we define 
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Figure 16. Contact number correlations for spheres in contact at (f>RCP- The error bars are standard 
deviations from the mean. The solid line are fits of Equation I I37I I to the data represented by closed symbols, 
open symbols are data omitted from fitting due to low occurrence. The dashed line is the prediction of an 
uncorrelated packing from Equation Il35|l . The data plotted in each panel is: (a) lognormal = 0.40; (b) 
Gaussian aj^ = 0.44; (c) uniform cr^ = 0.44; (d) bidisperse; (e) monodisperse. In (f) the fit parameter b is 
plotted as a function cr^. The data in (f) is labelled the same as in Figure[3la). 
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Znniz), the average contact number of particles that are in contact with a par- 
ticles that has z contacts. Znn{z) is analogous to the quantity m{n), the average 
number of neighbours neighbouring a cell with n neighbours, in the Aboav-Weaire 
law, which is well studied for two dimensional cellular structures [s^, ^391]. Recently 
it has been shown that two dimensional disc packings exhibit nearest neighbour 
correlations in the contact network [33] similar to the anti-correlations observed in 
cellular structures. 



In our previous work [33(] we proposed an analogue of the Aboav-Weaire law for 
cellular structures to describe these correlation in two dimensional disc packings. 
This analogue draws upon a counting argument [s^ that leads to a sum rule that 
is exact and independent of dimension, 

- z'^)P{z) = 0. (34) 



Znn is then a function of z that must satisfy Equation (j34p . For uncorrelated 
packings Znn{z) is a constant (= Znn) which is independent of z. Using Equation 
(j34p we can show that it is given by 

2 

Znn = {z) + ^. (35) 

Figure [TU] shows Z„„(z) for various size distributions. All distributions exhibit 
clear anti-correlations, namely that particles with few contacts are surrounded 
by particles with many contacts and vice versa. However, the deviations from the 
uncorrelated prediction Znn is usually less than 10%, which may explain the reason 
why the granocentric approach works well. 

While there is currently no theoretical prediction for these correlations, one can 
find an empirical equation for Znn{z) based on a series expansion of Znni^z) in 
terms of the moments of P{z). In order to ensure that the sum rule (Equation 
(134p ) is satisfied the series takes the form of 

{Znn - {Z))Z - al = {z' " (/)) , (36) 

i=l 

where the Cj's are arbitrary constants. If q = for z > 1, one recovers the Aboav- 
Weaire law for cellular structures. For 2D packings we found that the data was 



well described by only making the second term non-zero [33(] which leads to a one 
parameter fit: 

z 

where b = C2- This empirical equation agrees well with the correlations we find in 
2D [33] and 3D packings as shown in Figure [T6l 

The fit parameter b does not depend on the shape of the size distribution but 
only on the width a a as shown in Figure [T6lffl. Note that all size distributions 
regardless of shape or width exhibit these anti-correlations in the contact network. 



4.2. Nearest neighbour size correlations 



In the following we investigate spatial correlations between the size of particles in 
our polydisperse packings in two and three dimensions. Another way to look at this 
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Figure 17. (a) Correlations between size of spheres in contact. All symbols represent the same size distri- 
butions as in Figure FlBl except: ( ) Gaussian ua = 0.27. The solid lines are fits to the data using Equation 
II40I I with the dashed lines being the uncorrelated prediction calculated from Equation I I44I1 . The inset 
shows the fit parameter ui as a function of a a . The data in the inset is labelled the same as in FigurelSia). 
(b) Correlations between size of discs in contact in two dimensions rescaled by the predicted uncorrelated 
radius Run- Five different size distributions are plotted; (0) bidisperse, radius ratio 1:1.4; (□) uniform 
an = 0.17; ( ) Gaussian an = 0.14; ( i) lognormal an = 0.20; (v) lognormal an = 0.30; (>) lognormal 
an = 0.40 



question is as follows: Is the size distribution of particles neighbouring a central 
particle of radius rc just equal to the global size distribution P{r) as assumed in 
the granocentric approach? 

In order to explore potential size correlations, ^„„(a) is defined as the average 
normalised surface area of all particles in contact with a particle with area a. 
Figure [TTlf a) . shows Ann versus a for four different size distributions. This result 
clearly indicates spatial anti-correlations in the particle size and are not consistent 
with the uncorrelated prediction Ann which we discuss below. On average larger 
particles are surrounded by smaller particles and vice versa. 

The same counting argument that is used to formulate Equation (134p can be 
applied to the particle size. Analogous to Equation ([34ll . Ann must satisfy the 
following relation: 

/ Ann{a){z\a)P{a)da = / a{z\a)P{a)da. (38) 
Jo Jo 
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Using a series expansion in terms of the moments of the area distribution P{a) 
that satisfies Equation ([38|) we obtain 



. . ^ _ /o°° a{z\a)P{a)da + ^. u;i(a' - (g^)) 

{z\a} 

where Wi are arbitrary constants. Previously we have shown that {z\a) was weh 
described by a hnear relation (Equation ([9])). Substituting this expression into 
Equation (j39p and keeping only the first term in the expansion (i = 1) yields, 

{z)+jal + w{a-l) 

where u; = wi. This is a one parameter fit to the data, since the value of 7 is 
independent of polydispersity. The fits shown in Figure [TTlfa) agree well with the 
data for all the size distributions we considered. The inset shows the corresponding 
values of w which mostly depend on the width but not the shape of the size dis- 
tribution. Therefore, Equation ()40p provides a good description of the correlations 
between the size of particles in 3D disordered packings. 

In the absence of correlations A„„,(a) becomes a constant (= Ann) which we can 
derive from Equation (ISSp : 



00 



Ann= I a{ ^P{a) ] da. (41) 



V 

Substituting Equation ([9]) into the previous expression, it can be simplified to 

Ann = l + T^CTa, (42) 

where 7 and (z) are constant at 4>rcp- Therefore, Ann depends on aA only. This 
relation is plotted in Figure [18] and agrees well with the data. Wider distributions 
give rise to a larger Ann- 
Alternatively, the sum rule analogue for particle size can be written as the sum 

Y,{Ann\z)zP{z) = z{a\z)P{z), (43) 



where {Ann\z) is the average area of all particles in contact with a particle that 
has z contacts. In the absence of correlations {Ann\z) reduces to 

(^E^H^)^(^)- (44) 



This expression for Ann is equivalent to Equation (j4ip as can be seen from Equation 

m- 

We previously reported [s^ a sum rule for nearest neighbour size correlations in 
two dimensions 

oo POO 

Rnn{f){z\r)P{r)dr = I r{z\r)P{r)dr. (45) 
Jo 
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Figure 18. The uncorrelated prediction Ann for a range of different size distributions in three dimensions. 
The dashed line is the uncorrelated prediction of Ann calculated from Equation 1142 H . Inset: The uncorre- 
lated prediction Rnn for a range of different size distribution in two dimensions. The dashed line is the 
uncorrelated prediction of Knn calculated from Equation I I47II . The data is labelled the same as in Figure 
Ela). 



where Rnnif) is the average normahsed radius in contact with a particle of radius r. 
The data is shown in Figure [TTlfb). In contrast to the results for three dimensional 
sphere packings, the size distribution appears to be spatially uncorrelated. Rnn is 
well described by 



R„ 



°° ^ (z\r) 

r I VT^P(r) ) dr, 



for all polydipsersities. Using the the linear fit {z\r) this reduces to 

Rnn 



1 _L '^20 2 
1 + -T^(^R, 



(46) 



(47) 



where 72D and (z) are constant at (pRCP- 
Fi gure 1181 shows that both Ann and Rnn 

increase with a a and an, respectively. 
The data is well described by Equations ()32|) and ([I7|) . The reason for both Ann 
and Rnn to be greater than 1 is the fact larger particles have more contacts and 
therefore appear more frequently in the first neighbour shell. 

These results affect the central assumption in the granocentric approach, namely 
that the size distribution of the neighbour shell is equivalent to the global P(r). 
Even in the absence of correlations, as for 2D packings, the average particle radius 
of the first neighbour shell is larger than 1 - the mean radius of the global size 
distribution. When correlations are neglected, the effective size distribution P{rnn) 
of particles in the (contacting) neighbour shell is simply 



Pirnn) 



{z\r) 



P(r) 



(48) 
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which fohows from Equation (|46p . Consequently, the assumption that the size dis- 
tribution of the neighbouring particles is governed by P{r) becomes progressively 
worse for broader distribution. 

In practice, the difference between the weighted size distribution shown above 
and P{r) is small for moderate polydispersities and therefore does not measurably 
affect the outcome of the predictions made by the granocentric model. 



Conclusions 

We have shown that a surprising number of features of frictionless packings are 
insensitive to polydispersity. Our key result is the universal correlations we observe 
between size and contact number of a particle that are independent of the shape 
and width of the size distribution. This holds in both two and three dimensions 
and allows a mean field formulation of the granocentric model. The contact number 
distributions emerging from the model agree well with the data for a wide range 
of polydispersities. The two parameters that appear in the model do not depend 
on polydispersity either. 

In passing we note that the random close packing density excluding the rattlers 
remains unchanged for a wide range of discrete and continuous size distributions in 
both dimensions. This holds even for packings that contain up to 30% of rattlers. 

Despite the success of the granocentric approach, we note that a central assump- 
tion in this model does not hold. We find that all packings exhibit anti-correlations 
in the contact network. In addition, the particle sizes are anti-correlated, but only 
in 3D packings. For three dimensional packings we can therefore conclude that on 
average smaller particles which have less contacts are surrounded by larger particles 
that have more contacts. 

Nevertheless, the granocentric model, while not exact, yields good predictions 
since these correlations are typically weak. 
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